package scu.maqiang.possion;

import java.util.Arrays;

import scu.maqiang.fes.BCType;
import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2Q81;
import scu.maqiang.mesh.Mesh2Q8;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class Possion2DQ81 {

	public static void main(String[] args) {
        int N = 8, label = 10;
        double ep = 1.0 / N;
        
        //ʹ��lambda���ʽ�����ί��
        ScalarFunc regionFunc = (xy, llabel, param) -> {
            double xCell = xy[0] * N - (int)(xy[0] * N);
            double yCell = xy[1] * N - (int)(xy[1] * N);
            return (Math.abs(xCell - 0.5) < 0.25 && Math.abs(yCell - 0.5) < 0.25)? 1.0:-1.0;
        };

		Mesh2Q8 mesh2D = new Mesh2Q8().square(128, 128);
        mesh2D.setDomainLabel(regionFunc, null, label);
        
        FES2Q81 fs = new FES2Q81(mesh2D);
        SRMatrix A = new SRMatrix(fs.GetNdof());

        ScalarFunc regionCoef = (xyz, relabel, param) -> relabel == label ? 0.01 : 1.0;
        ScalarFunc ConstCoef = (xyz, relabel, param) -> 1.0e0;
        ScalarFunc[] gcf = {regionCoef};
        gcf[0] = ConstCoef;
        
        ScalarFunc coefFunc = (xyz, relabel, param) -> 2 * xyz[0] * (1 - xyz[0]) + 2 * xyz[1] * (1 - xyz[1]);
        fs.assembleStiff(gcf, null, BVPType.COMMON, A);
        //fs.assembleHeatStiff(new double[] {1.0}, BVPType.COMMON, A);
        //A.sort();
        //System.out.println(A);
//        System.out.println(A.makeEntrySymmetric());
        double[] RHS = new double[fs.GetNdof()];
        //GeneralCoefFunc ConstCoef = (xyz, relabel, param) -> 1.0;
//        coefFunc = (xy, relabel, param) ->  {
//        	double a = 400.0;
//        	double r = Math.sqrt((xy[0] - 0.5) * (xy[0] - 0.5) + (xy[1] - 0.50) * (xy[1] - 0.5));
//        	return 4 * 400 * 400 * (1 - a * r * r) * Math.exp(-400.0 * r * r);
//        };
        gcf[0] = coefFunc;
        fs.assembleSource(gcf, null, BVPType.COMMON, RHS);
//        System.out.println(MVO.toString(RHS));
        A.makeEntrySymmetric();
        fs.applyBCZero(A, RHS, Direct.X, BCType.MBN, 1, 2, 3, 4);
        
        //System.out.println(A);
        //System.out.println(Arrays.toString(RHS));
        
        

        double[] x = new double[fs.GetNdof()];
        
        //IterSSolver solver = new IterSSolver(1.0e-15, 1000000);
        
        //CSRMatrix mat = new CSRMatrix(A);
        //solver.PCGSSOR(A, RHS, x, 1.5e0, 1);
        //Arrays.fill(x, 0.0);
        //solver.PCGSSOR_2(A, RHS, x, 1.5e0, 1);
        
        //System.out.println(A);
        
        //System.out.println(mat);
        
        //Arrays.fill(x, 0.0);
        //solver.CG(mat, RHS, x, 1);
        
        //Arrays.fill(x, 0.0);
        //solver.CG(A, RHS, x, 1);
        
        Arrays.fill(x, 0.0);
        NewIterSSolver newSolver = new NewIterSSolver(A);
        newSolver.PCGSSOR(RHS, x, 1.5e0, 1);

        mesh2D.toTecplot("Solutions.dat", x);
 
        double[][] gradX = new double[2][fs.GetNdof()];
        fs.computeGradient2(x, gradX);
        mesh2D.toTecplot("Solutions.dat", x);
        mesh2D.toTecplot("SolutionsGradient.dat", gradX);
        System.out.println("Gradient: ");
        System.out.println("Dx min: " + MVO.min(gradX[0]) + "\tmax: " + MVO.max(gradX[0]));
        System.out.println("Dy min: " + MVO.min(gradX[1]) + "\tmax: " + MVO.max(gradX[1]));
        
        double[][] sdX = new double[3][fs.GetNdof()];
        fs.computeSecondDerivative(x, sdX);
        System.out.println("Second Derivative: ");
        System.out.println("Dxx min: " + MVO.min(sdX[0]) + "\tmax: " + MVO.max(sdX[0]));
        System.out.println("Dxy min: " + MVO.min(sdX[1]) + "\tmax: " + MVO.max(sdX[1]));
        System.out.println("Dyy min: " + MVO.min(sdX[2]) + "\tmax: " + MVO.max(sdX[2]));
        mesh2D.toTecplot("SolutionsSD.dat", sdX);
        
        fs.computeFirstAndSecondDerivative(x, gradX, sdX);
	}

}
